#include <iostream.h>
#include <fstream.h>
#include <stdio.h>
#include <math.h>
#include <datafile.h>
#include <cmlinset.h>
#include <spctmdp.h>
#include <ctmdpvi.h>
#include <ctvianu.h>
#include <ctviau.h>
#include <ctviasu.h>
#include <ctstclas.h>

#include <ctmcemb.h>
#include <cmstprob.h>
#include <timer.h>
#include <vector>

//const char *PATH = "c:\\solon\\cpp\\OpticalWDM2\\";
// PATH para as bibliotecas do Modesto(Biblioteca implementada pelo Prof. Solon Ven�ncio);
const char *PATH = "D:\\P�s_Gradua��o\\Disserta��o\\Implementa��o\\OpticalWDM - Rede Heterog�nea - Vers�o 0.4\\OpticalWDM2\\";

// Declara��o de Vari�veis;

const int N1 = 0, N2 = 1, N3 = 2, EV = 3, // componentes do sistema;
        S1 = 0, S2 = 1, S3 = 2, A1 = 3, A2 = 4, A3 = 5, // eventos;
        NOA = 0, REJ = 0, ACC = 1; // a��es;
//N1 = chamadas da classe 1, N2 = chamadas da classe 2, N3 = chamadas da classe 3, EV = Evento;
//S1 = T�rmino uma chamada da classe 1, S2 = T�rmino uma chamada da classe 2, S3 = T�rmino uma chamada da classe 3,
//A1 = Chegada uma chamada da classe 1, A2 = Chegada uma chamada da classe 2, A3 = Chegada uma chamada da classe 3,
//NOA = Nenhuma A��o, REJ = Rejeita da Chamada, ACC = Aceita uma Chamada;

UNSIGNED nWL; //considera os valores positivos - 0 a 255;
FLOAT lambCon, muCon, lambClass1, lambClass2, lambClass3,
lambInNextNodeClass1, lambInNextNodeClass2, lambInNextNodeClass3,
reward1, reward2, reward3, UBReward;
//lambCon: taxa de chegada das chamadas em cada conex�o : 6;
//muCon: taxa de processamento em cada conex�o : 1;
//lambClass1: taxa do fluxo que chega no n� de observa��o;
//lambClass2: taxa do fluxo que passa pelo n� de observa��o;
//lambClass3: taxa do fluxo que sai do n� de observa��o;
//lambInNextNodeClass1: taxa do fluxo que chega no pr�ximo n�;
//lambInNextNodeClass2: taxa do fluxo que passa pelo pr�ximo n�;
//lambInNextNodeClass3: taxa do fluxo que sai do pr�ximo n�;
//reward1: retorno para a Classe 1 (chamadas que chegam no n� de observa��o);
//reward2: retorno para a Classe 2 (chamadas que passam pelo n� de observa��o);
//reward3: retorno para a classe 3 (chamadas que partem do n� de observa��o);
//UBReward: par�metro a ser descoberto;

//Declara��o de vari�veis
//FLOAT  probLoss1=0.0, probLoss2=0.0, probLoss3=0.0,
FLOAT probBlock1, probBlock2, probBlock3;

//FLOAT newProbLoss1,  newProbLoss2,  newProbLoss3;
FLOAT IncomingNetRateInNode, OutgoingNetRateInNode,
CrossingInNodeRate, CrossingOutNodeRate;

/*  Array2<FLOAT> lambIn(4, 4, 0.0),
                lambOut(4, 4, 0.0), // D�vida!!!!
                lambInInNextNode(4, 4, 0.0);
 */
int numNos = 4; //numero de n�s

Array2<FLOAT> lambIn(numNos, numNos, 0.0),
lambOut(numNos, numNos, 0.0),
lambInInNextNode(numNos, numNos, 0.0),
ProbBlock(numNos, numNos, 0.0);

//lambIn: Taxa de chegada no n�;
//lambOut: Taxa de sa�da no n�;
//lambInInNextNode: Taxa de entrada no pr�ximo n�;



UNSIGNED maxIter, optimizationMethod, completeOutput, useLastIterationData, //considera os valores positivos - 0 a 255;
objFuncType;
//maxIter: n�mero m�ximo de itera��es;
//optimizationMethod: determina o m�todo de otimiza��o;
//completeOutput: sa�da completa;
//useLastIterationData: utiliza os dados da �ltima itera��o;
//objFuncType: determina o tipo da fun��o objetivo;

FLOAT tauFactor, unifFactor, //pormenores matem�ticos;
limProbRelaxFactor, rateUpdateRelaxFactor,
viPrecision, lpPrecision;
//tauFactor: VI Tau factor (Tau factor) : 1.0;
//unifFactor: VI uniformization factor : 0.1;
//limProbRelaxFactor: limit Probability Relaxation Factor: 1    0.5;
//rateUpdateRelaxFactor: rates Between nodes Relaxation Factor : 1.0;
//viPrecision: Precision for Value Iterator Algorithm : 1.0e-8;
//lpPrecision: Precision for Limit Probability Algorithm : 1.0e-10;
//Onde VI = Algoritmo de Itera��o de Valores;

Array<FLOAT> Val, limitProb;
//Val: vetor tipo float para armazenar valores; (muito obscuro ainda!!)
//limitProb: vetor tipo float para armazenar as probabilidades limites;

UNSIGNED nReservedWL1, nReservedWL2, nReservedWL3; //considera os valores positivos - 0 a 255;
//nReservedWL1: n�mero de comprimentos de onda reservados para as chamadas que chegam ao n� de observa��o;
//nReservedWL2: n�mero de comprimentos de onda reservados para as chamadas que passam pelo n� de observa��o;
//nReservedWL3: n�mero de comprimentos de onda reservados para as chamadas que partem do n� de observa��o;

vector<vector<double> > Matrix_lambCon(numNos, vector<double>(numNos));
vector<vector<vector<double> > > Set_ProbBlock(numNos, Matrix_lambCon);



//Classe OpticalWDMCondition
//Heran�a entre classess. A classe OpticalWDMCondition herda as caracter�sticas
//da classe Condition;

class OpticalWDMCondition : public Condition {
protected:
    UNSIGNED W; //N�mero de comprimentos de onda;
public:

    OpticalWDMCondition(UNSIGNED nWL) {
        W = nWL;
    }

    OpticalWDMCondition(OpticalWDMCondition & OC) {
        W = OC.W;
    }

    virtual Condition * getACopy() {
        return new OpticalWDMCondition(*this);
    }
    virtual int isVerifiedFor(MultiLinearSet::Element & e);

    virtual void printOn(ostream &os) {
        os << "Optical WDM condition";
    }

    virtual void printOnIdented(ostream &os, int ident) {
        identOStream(os, ident);
        printOn(os);
    }

    virtual setClass isA() {
        return auxiliar0CondClass;
    }

    virtual void clean() {
    }

    virtual ~OpticalWDMCondition() {
    } //desconstrutor;
};
//m�todo                                                                //OBJETO DA CLASSE ELEMENT, DA CLASSE MULTILINEARSET

int OpticalWDMCondition::isVerifiedFor(MultiLinearSet::Element & se) {
    UNSIGNED ev = se.getCoordValue(EV), //Evento;
            n1 = se.getCoordValue(N1), //Chamadas que chegam no n� de observa��o;
            n2 = se.getCoordValue(N2), //Chamadas que passam pelo n� de observa��o;
            n3 = se.getCoordValue(N3); //Chamadas que partem do n� em dire��o a outro;
    if (ev == S1 && n1 == 0 || //
            ev == S2 && n2 == 0 ||
            ev == S3 && n3 == 0 ||
            n1 + n2 > W ||
            n3 + n2 > W)
        return 0;
    return 1;
}
//objeto = endere�o;

int getNActions(MultiLinearSet::Iterator & si) {
    UNSIGNED ev = si().getCoordValue(EV);
    if (ev == S1 || ev == S2 || ev == S3)
        return 1;
    else {
        UNSIGNED n1 = si().getCoordValue(N1),
                n2 = si().getCoordValue(N2),
                n3 = si().getCoordValue(N3);
        if (ev == A1 && n1 + n2 == nWL ||
                ev == A2 && (n1 + n2 == nWL || n3 + n2 == nWL) ||
                ev == A3 && n3 + n2 == nWL)
            return 1;
    }
    return 2;
}
//m�todo retorna um objeto do tipo index;

INDEX maxAcceptanceLevel(MultiLinearSet::Iterator & si, INDEX nwlTotal, INDEX c) {
    UNSIGNED n1 = si().getCoordValue(N1),
            n2 = si().getCoordValue(N2),
            n3 = si().getCoordValue(N3);
    switch (c) {
        case 1: return nwlTotal - n2;
        case 2: return nwlTotal - (n1 > n3) ? n1 : n3;
        case 3: return nwlTotal - n2;
    }
    return 0;
}

INDEX maxAcceptanceLevel(INDEX nwlTotal, INDEX n1, INDEX n2, INDEX n3, INDEX c) {
    switch (c) {
        case 1: return nwlTotal - n2;
        case 2: return nwlTotal - (n1 > n3) ? n1 : n3;
        case 3: return nwlTotal - n2;
    }
    return 0;
}

//m�todo;                //ponteiro - refer�ncia

int getDataFromFile(char * fName) {
    DataFile File(fName); // objeto
    char ch;
    if (!File.isOk()) {
        cerr << "\nCan't not open " << fName;
        cin >> ch;
        exit(1);
    }
    cout << "\nReading file " << fName;

    File.readMsg(nWL, "number of available wavelengths");
    File.readMsg(lambCon, "call on rate in each conexion");
    File.readMsg(muCon, "call off rate in a each conexion");

    File.readMsg(objFuncType, "objective function type");
    File.readMsg(reward1, "reward rate of a class 1 user");
    File.readMsg(reward2, "reward rate of a class 2 user");
    File.readMsg(reward3, "reward rate of a class 3 user");
    if (objFuncType == 3 || objFuncType == 4) {
        reward1 = 1.0;
        reward2 = 1.0;
        reward3 = 1.0;
    }
    UBReward = (reward1 > reward2) ? reward1 : reward2;
    UBReward = (UBReward > reward3) ? UBReward : reward3;
    UBReward = nWL*UBReward;
    File.readMsg(optimizationMethod, "optimization method");
    File.readMsg(nReservedWL1, "Given number of wavelenght for class 1");
    File.readMsg(nReservedWL2, "Given number of wavelenght for class 2");
    File.readMsg(nReservedWL3, "Given number of wavelenght for class 3");
    File.readMsg(viPrecision, "Precision for Value Iterator Algorithm");
    File.readMsg(lpPrecision, "Precision for Limit Probability Algorithm");
    File.readMsg(limProbRelaxFactor, "Limit Probability Relaxation Factor");
    File.readMsg(tauFactor, "Tau factor");
    File.readMsg(unifFactor, "Uniformization factor");
    File.readMsg(rateUpdateRelaxFactor, "Rates Between nodes Relaxation Factor");
    File.readMsg(maxIter, "Max number of iterations for VI algorithm");
    File.readMsg(useLastIterationData, "Use Last Iteration Data Flag");
    File.readMsg(completeOutput, "Complete output flag");
    return 1;
}

void showDataOn(ostream & os) {
    os << "\n\n******************** Optical WDM Model ***********************\n\n";
    os << "\n            Optimal Resource Allocation and Fairness";
    os << "\n              Control in All-Optical WDM Networks";
    os << "\n               K.Mosharaf, J.Talim, I.Lambadaris";
    os << "\nIEEE JOURNAL ON SELECTED AREAS IN COMMUNICATIONS, V23, N8, 2005\n\n";


    os << "\n\nnumber of available wavelengths: " << nWL;
    os << "\n\nCall on rate in each conexion:   " << lambCon;
    os << "\nCall off rate in each conexion: " << muCon;


    os << "\nobjective function type:"
            << "\n   Maximize: ";
    switch (objFuncType) {
        case 1: os << "node weighted reward rate of each two-nodes wavelenght used by each class user";
            break;
        case 2: os << "node weighted throughput each class user";
            break;
        case 3: os << "network reward rate of each two-nodes wavelenght";
            break;
        case 4: os << "VI with Total Information and semi-uniformization";
            break;
        default: cerr << "\nnetwork throughput: " << objFuncType;
    }
    if (objFuncType == 1 || objFuncType == 2) {
        os << "\n\nreward rate of a class 1 user:   " << reward1;
        os << "  \nreward rate of a class 2 user:   " << reward2;
        os << "  \nreward rate of a class 3 user:   " << reward3;
        os << "  \nUpper Bound for reward rate of the system: " << UBReward;
    }
    os << "\n\nOptimization method: ";
    switch (optimizationMethod) {
        case 1: os << "VI with Total Information and uniformization";
            break;
        case 2: os << "Given number of wavelenght for each class in node";
            break;
        case 3: os << "VI with Total Information and no uniformization";
            break;
        case 4: os << "VI with Total Information and semi-uniformization";
            break;
        default: cerr << "\nInvalid value for information type: "
                    << optimizationMethod;
    }
    if (optimizationMethod == 2) {
        os << "\n   Class 1 (quiting the ring in the node):  " << nReservedWL1;
        os << "\n   Class 2 (crossing the node):             " << nReservedWL2;
        os << "\n   Class 3 (entering the ring in the node): " << nReservedWL3;
    }
    os << "\n\nPrecision for Value Iterator Algorthm:   " << viPrecision;
    os << "\nPrecision for Limit Probability Algorithm: " << lpPrecision;
    os << "\nLimit Probability  Relaxation Factor:      " << limProbRelaxFactor;
    os << "\nTau factor:                                " << tauFactor;
    os << "\nUniformization factor:                     " << unifFactor;
    os << "\nRates Between nodes Relaxation Factor:     " << rateUpdateRelaxFactor;
    os << "\nMax number of iterations f or VI algorithm : " << maxIter;
    os << "\nUse Last Iteration Data Flag (0, 1)  : " << useLastIterationData;
    os << "\nComplete output flag (0, 1)          : " << completeOutput;
    os << "\n\n";
    os.flush();
}

void printState(ostream & os, MultiLinearSet::Element & se) {
    os << "(" << setw(2) << se.getCoordValue(N1) << ", "
            << setw(2) << se.getCoordValue(N2) << ", "
            << setw(2) << se.getCoordValue(N3) << ", ";
    switch (se.getCoordValue(EV)) {
        case S1: os << "S1)";
            break;
        case S2: os << "S2)";
            break;
        case S3: os << "S3)";
            break;
        case A1: os << "A1)";
            break;
        case A2: os << "A2)";
            break;
        case A3: os << "A3)";
            break;
        default: os << "??";
    }
}

void printAction(ostream & os, MultiLinearSet::Element & se, int action) {
    UNSIGNED ev = se.getCoordValue(EV);
    if (action == NOA) {
        if (ev == S1 || ev == S2 || ev == S3)
            os << "NOA";
        else {
            UNSIGNED n1 = se.getCoordValue(N1),
                    n2 = se.getCoordValue(N2),
                    n3 = se.getCoordValue(N3);
            if (ev == A1 && n1 + n2 == nWL ||
                    ev == A2 && (n1 + n2 == nWL || n3 + n2 == nWL) ||
                    ev == A3 && n3 + n2 == nWL)
                os << "rej";
            else
                os << "REJ";
        }
    } else
        os << "ACC";


}

void printTransition(ostream & os, MultiLinearSet::Element & sf,
        MultiLinearSet::Element & st, int a, char * rate) {
    printState(os, sf);
    os << "  ";
    printAction(os, sf, a);
    os << "  ";
    os.flush();
    printState(os, st);
    os << "  " << rate << "\n";
}

void reactToAction(MultiLinearSet::Element & sf, int a,
        MultiLinearSet::Element & sr) {
    sr = sf;
    switch (sf.getCoordValue(EV)) {
        case S1: sr.decCoordValue(N1);
            break;
        case S2: sr.decCoordValue(N2);
            break;
        case S3: sr.decCoordValue(N3);
            break;
        case A1: if (a == ACC) sr.incCoordValue(N1);
            break;
        case A2: if (a == ACC) sr.incCoordValue(N2);
            break;
        case A3: if (a == ACC) sr.incCoordValue(N3);
            break;
    }
}

void setFixedReservationPolicy(ConstrMultiLinearSet & CS, Policy & pol,
        INDEX nwlTotal, INDEX nwl1, INDEX nwl2, INDEX nwl3) {
    UNSIGNED ev, n1, n2, n3, nwl1or3;
    MultiLinearSet::Iterator si(CS);
    nwl1or3 = (nwl1 < nwl3) ? nwl3 : nwl1;
    for (si.reset(); !si.isDone(); si.next()) {
        ev = si().getCoordValue(EV);
        if (ev == S1 || ev == S2 || ev == S3)
            pol.setAction(si, NOA);
        else {
            n1 = si().getCoordValue(N1);
            n2 = si().getCoordValue(N2);
            n3 = si().getCoordValue(N3);
            if (ev == A1 && (n1 + n2 == nWL || n1 >= nwlTotal - nwl2) ||
                    ev == A2 && (n1 + n2 == nWL || n3 + n2 == nWL || n2 >= nwlTotal - nwl1or3) ||
                    ev == A3 && (n3 + n2 == nWL || n3 >= nwlTotal - nwl2))
                pol.setAction(si, REJ);
            else
                pol.setAction(si, ACC);
        }
    }
}

INDEX getFixedReservationPolicyAction(MultiLinearSet::Iterator & si,
        INDEX nwlTotal, INDEX nwl1, INDEX nwl2, INDEX nwl3) {
    UNSIGNED ev;
    ev = si().getCoordValue(EV);
    if (ev == S1 || ev == S2 || ev == S3)
        return NOA;
    UNSIGNED n1, n2, n3, nwl1or3;
    nwl1or3 = (nwl1 < nwl3) ? nwl3 : nwl1;
    n1 = si().getCoordValue(N1);
    n2 = si().getCoordValue(N2);
    n3 = si().getCoordValue(N3);
    if (ev == A1 && n1 >= nwlTotal - nwl2 ||
            ev == A2 && n2 >= nwlTotal - nwl1or3 ||
            ev == A3 && n3 >= nwlTotal - nwl2)
        return REJ;
    return ACC;
}

void calcMinReservationLevelsForAPolicy(ConstrMultiLinearSet & CS, Policy & pol,
        INDEX nwlTotal, INDEX & minNWL1, INDEX & minNWL2, INDEX & minNWL3) {
    UNSIGNED ev, n1, n2, n3;

    minNWL1 = minNWL3 = minNWL2 = nwlTotal;
    MultiLinearSet::Iterator si(CS);
    for (si.reset(); !si.isDone(); si.next()) {
        ev = si().getCoordValue(EV);
        if (ev == A1 || ev == A2 || ev == A3) {
            n1 = si().getCoordValue(N1);
            n2 = si().getCoordValue(N2);
            n3 = si().getCoordValue(N3);
            if (ev == A1 && pol.getAction(si) == ACC && minNWL2 >= nwlTotal - n1 - 1)
                minNWL2 = nwlTotal - n1 - 1;
            if (ev == A2 && pol.getAction(si) == ACC) {
                if (minNWL1 >= nwlTotal - n2 - 1)
                    minNWL1 = nwlTotal - n2 - 1;
                if (minNWL3 >= nwlTotal - n2 - 1)
                    minNWL3 = nwlTotal - n2 - 1;
            }
            if (ev == A3 && pol.getAction(si) && minNWL2 >= nwlTotal - n3 - 1)
                minNWL2 = nwlTotal - n3 - 1;
        }
    }
}

void PolicyAnalysis(ostream & os, ConstrMultiLinearSet & CS,
        Policy & pol, INDEX nwlTotal) {
    UNSIGNED i, j, a, ev, n1, n2, n3, ni, no;
    const UNSIGNED ACCEPT1 = 0x1, REJECT1 = 0x2,
            ACCEPT2 = 0x4, REJECT2 = 0x8,
            ACCEPT3 = 0x10, REJECT3 = 0x20;
    bool ok;
    Array2<int> LBRej1(nwlTotal + 1, nwlTotal + 1, nwlTotal + 1),
            LBRej2(nwlTotal + 1, nwlTotal + 1, nwlTotal + 1),
            LBRej3(nwlTotal + 1, nwlTotal + 1, nwlTotal + 1),
            UBAcc1(nwlTotal + 1, nwlTotal + 1, -1),
            UBAcc2(nwlTotal + 1, nwlTotal + 1, -1),
            UBAcc3(nwlTotal + 1, nwlTotal + 1, -1);

    Array2<UNSIGNED> InOutAction(nwlTotal + 1, nwlTotal + 1, 0x0);

    MultiLinearSet::Iterator si(CS);
    for (si.reset(); !si.isDone(); si.next()) {
        ev = si().getCoordValue(EV);
        n1 = si().getCoordValue(N1);
        n2 = si().getCoordValue(N2);
        n3 = si().getCoordValue(N3);
        ni = n1 + n2;
        no = n2 + n3;
        a = pol.getAction(si);
        switch (ev) {
            case A1:
                if (a == ACC) {
                    if (UBAcc1(n2, n3)< (int) n1)
                        UBAcc1(n2, n3) = n1;
                    InOutAction(ni, no) = InOutAction(ni, no) | ACCEPT1;
                }
                if (a == REJ) {
                    if (LBRej1(n2, n3)> (int) n1)
                        LBRej1(n2, n3) = n1;
                    InOutAction(ni, no) = InOutAction(ni, no) | REJECT1;
                }
                break;
            case A2:
                if (a == ACC) {
                    if (UBAcc2(n1, n3)< (int) n2)
                        UBAcc2(n1, n3) = n2;
                    InOutAction(ni, no) = InOutAction(ni, no) | ACCEPT2;
                }
                if (a == REJ) {
                    if (LBRej2(n1, n3)> (int) n2)
                        LBRej2(n1, n3) = n2;
                    InOutAction(ni, no) = InOutAction(ni, no) | REJECT2;
                }
                break;
            case A3:
                if (a == ACC) {
                    if (UBAcc3(n1, n2)< (int) n3)
                        UBAcc3(n1, n2) = n3;
                    InOutAction(ni, no) = InOutAction(ni, no) | ACCEPT3;
                }
                if (a == REJ) {
                    if (LBRej3(n1, n2)> (int) n3)
                        LBRej3(n1, n2) = n3;
                    InOutAction(ni, no) = InOutAction(ni, no) | REJECT3;
                }
                break;
        }
    }
    os << "\n\n\nMax acceptance level for Class 1 arrivals:";
    os << "\n\nn2\\n3";
    for (j = 0; j <= nwlTotal; ++j)
        os << "  " << setw(2) << j;
    os << "\n";
    for (i = 0; i <= nwlTotal; ++i) {
        os << "\n" << setw(2) << i << ":  ";
        for (j = 0; j <= nwlTotal; ++j) {
            os << "  ";
            if (UBAcc1(i, j) >= 0)
                os << setw(2) << UBAcc1(i, j);
            else
                os << "--";
        }
    }
    os << "\n\nMin rejection level = Max acceptance level+1 for Class 1 arrivals";
    for (i = 0, ok = false; i <= nwlTotal; ++i) {
        for (j = 0; j <= nwlTotal; ++j) {
            if (!((UBAcc1(i, j) + 1 == LBRej1(i, j)) ||
                    (UBAcc1(i, j) == -1 && LBRej1(i, j) == nwlTotal + 1))) {
                if (!ok) {
                    os << "\nEXCEPT FOR: ";
                    ok = true;
                }
                os << "\n    n2 = " << i << " e n3 = " << j;
            }
        }
    }

    os << "\n\n\nMax acceptance level for Class 2 arrivals:";
    os << "\n\nn1\\n3";
    for (j = 0; j <= nwlTotal; ++j)
        os << "  " << setw(2) << j;
    os << "\n";
    for (i = 0; i <= nwlTotal; ++i) {
        os << "\n" << setw(2) << i << ":  ";
        for (j = 0; j <= nwlTotal; ++j) {
            os << "  ";
            if (UBAcc2(i, j) >= 0)
                os << setw(2) << UBAcc2(i, j);
            else
                os << "--";
        }
    }
    os << "\n\nMin rejection level = Max acceptance level+1 for Class 2 arrivals";
    for (i = 0, ok = false; i <= nwlTotal; ++i) {
        for (j = 0; j <= nwlTotal; ++j) {
            if (!((UBAcc2(i, j) + 1 == LBRej2(i, j)) ||
                    (UBAcc2(i, j) == -1 && LBRej2(i, j) == nwlTotal + 1))) {
                if (!ok) {
                    os << "\nEXCEPT FOR: ";
                    ok = true;
                }
                os << "\n    n1 = " << i << " e n3 = " << j;
            }
        }
    }

    os << "\n\n\nMax acceptance level for Class 3 arrivals:";
    os << "\n\nn1\\n2";
    for (j = 0; j <= nwlTotal; ++j)
        os << "  " << setw(2) << j;
    os << "\n";
    for (i = 0; i <= nwlTotal; ++i) {
        os << "\n" << setw(2) << i << ":  ";
        for (j = 0; j <= nwlTotal; ++j) {
            os << "  ";
            if (UBAcc3(i, j) >= 0)
                os << setw(2) << UBAcc3(i, j);
            else
                os << "--";
        }
    }
    os << "\n\nMin rejection level = Max acceptance level+1 for Class 3 arrivals";
    for (i = 0, ok = false; i <= nwlTotal; ++i) {
        for (j = 0; j <= nwlTotal; ++j) {
            if (!((UBAcc3(i, j) + 1 == LBRej3(i, j)) ||
                    (UBAcc3(i, j) == -1 && LBRej3(i, j) == nwlTotal + 1))) {
                if (!ok) {
                    os << "\nEXCEPT FOR: ";
                    ok = true;
                }
                os << "\n    n1 = " << i << " e n2 = " << j;
            }
        }
    }
    os << "\n\n\nAcceptance for Class 1 arrivals in function of ni=n1+n2 " <<
            "and no=n2+n3:";
    os << "\n\nni\\no";
    for (j = 0; j <= nwlTotal; ++j)
        os << " " << setw(2) << j;
    os << "\n";
    for (i = 0; i <= nwlTotal; ++i) {
        os << "\n" << setw(2) << i << ":  ";
        for (j = 0; j <= nwlTotal; ++j) {
            os << "  ";
            if ((InOutAction(i, j) & ACCEPT1) && !(InOutAction(i, j) & REJECT1))
                os << "A";
            if ((InOutAction(i, j) & ACCEPT1) && (InOutAction(i, j) & REJECT1))
                os << "a";
            if (!(InOutAction(i, j) & ACCEPT1) && (InOutAction(i, j) & REJECT1))
                os << "R";
            if (!(InOutAction(i, j) & ACCEPT1) && !(InOutAction(i, j) & REJECT1))
                os << "-";
        }
    }
    os << "\n\n\nAcceptance for Class 2 arrivals in function of ni=n1+n2 " <<
            "and no=n2+n3:";
    os << "\n\nni\\no";
    for (j = 0; j <= nwlTotal; ++j)
        os << " " << setw(2) << j;
    os << "\n";
    for (i = 0; i <= nwlTotal; ++i) {
        os << "\n" << setw(2) << i << ":  ";
        for (j = 0; j <= nwlTotal; ++j) {
            os << "  ";
            if ((InOutAction(i, j) & ACCEPT2) && !(InOutAction(i, j) & REJECT2))
                os << "A";
            if ((InOutAction(i, j) & ACCEPT2) && (InOutAction(i, j) & REJECT2))
                os << "a";
            if (!(InOutAction(i, j) & ACCEPT2) && (InOutAction(i, j) & REJECT2))
                os << "R";
            if (!(InOutAction(i, j) & ACCEPT2) && !(InOutAction(i, j) & REJECT2))
                os << "-";
        }
    }
    os << "\n\n\nAcceptance for Class 3 arrivals in function of ni=n1+n2 " <<
            "and no=n2+n3:";
    os << "\n\nni\\no";
    for (j = 0; j <= nwlTotal; ++j)
        os << " " << setw(2) << j;
    os << "\n";
    for (i = 0; i <= nwlTotal; ++i) {
        os << "\n" << setw(2) << i << ":  ";
        for (j = 0; j <= nwlTotal; ++j) {
            os << "  ";
            if ((InOutAction(i, j) & ACCEPT3) && !(InOutAction(i, j) & REJECT3))
                os << "A";
            if ((InOutAction(i, j) & ACCEPT3) && (InOutAction(i, j) & REJECT3))
                os << "a";
            if (!(InOutAction(i, j) & ACCEPT3) && (InOutAction(i, j) & REJECT3))
                os << "R";
            if (!(InOutAction(i, j) & ACCEPT3) && !(InOutAction(i, j) & REJECT3))
                os << "-";
        }
    }

}

void go(ConstrMultiLinearSet & CS, ofstream & fileOut, FLOAT lambClass1,
        FLOAT lambClass2,
        FLOAT lambClass3,
        FLOAT & probBlock1,
        FLOAT & probBlock2,
        FLOAT & probBlock3) {
    char ch;
    Array<UNSIGNED> nAct(CS.getNElements());
    MultiLinearSet::Iterator sif(CS), sit(CS);

    for (sif.reset(); !sif.isDone(); sif.next())
        nAct[sif.getPosition()] = getNActions(sif);
    int i, j, n, a;
    if (completeOutput) {
        fileOut << "\n\nPossible States and actions:\n\n";
        for (sif.reset(); !sif.isDone(); sif.next()) {
            fileOut << setw(2) << sif.getPosition() << "   ";
            printState(fileOut, sif());
            for (a = 0; a < nAct[sif.getPosition()]; ++a) {
                fileOut << "   ";
                printAction(fileOut, sif(), a);
            }
            fileOut << "\n";
        }
    }
    Policy policy(CS);
    SpCTMkvDecProc Proc(CS, nAct);
    MultiLinearSet::Element elemAux(CS.getDimension());

    cout << "\n Generating transitions... \n";
    if (completeOutput)
        fileOut << "\n\nPOSSIBLE TRANSITIONS:\n\n";

    char ratestr[50];
    UNSIGNED n1, n2, n3, ev, n1t, n2t, n3t;
    for (sif.reset(); !sif.isDone(); sif.next()) {
        //n1 = sif().getCoordValue(N1);
        //n2 = sif().getCoordValue(N2);
        //n3 = sif().getCoordValue(N3);
        ev = sif().getCoordValue(EV);
        for (a = 0; a < Proc.getNActions(sif); ++a) {
            // elemAux = sif();
            // if(ev==S1)           elemAux.decCoordValue(N1);
            // if(ev==S2)           elemAux.decCoordValue(N2);
            // if(ev==S3)           elemAux.decCoordValue(N3);
            // if(ev==A1 && a==ACC) elemAux.incCoordValue(N1);
            // if(ev==A2 && a==ACC) elemAux.incCoordValue(N2);
            // if(ev==A3 && a==ACC) elemAux.incCoordValue(N3);
            reactToAction(sif(), a, elemAux);
            n1t = elemAux.getCoordValue(N1);
            n2t = elemAux.getCoordValue(N2);
            n3t = elemAux.getCoordValue(N3);
            if (n1t > 0) {
                elemAux.setCoordValue(EV, S1);
                sprintf(ratestr, "%d*mu", n1t);
                if (completeOutput) {
                    printTransition(fileOut, sif(), elemAux, a, ratestr);
                }
                if (!sit.goToElement(elemAux))
                    fileOut << "\n Error in transition";
                fileOut.flush();
                Proc.addRate(sif, sit, a, n1t * muCon);
            }
            if (n2t > 0) {
                elemAux.setCoordValue(EV, S2);
                sprintf(ratestr, "%d*mu", n2t);
                if (completeOutput) {
                    printTransition(fileOut, sif(), elemAux, a, ratestr);
                }
                if (!sit.goToElement(elemAux))
                    fileOut << "\n Error in transition";
                fileOut.flush();
                Proc.addRate(sif, sit, a, n2t * muCon);
            }
            if (n3t > 0) {
                elemAux.setCoordValue(EV, S3);
                sprintf(ratestr, "%d*mu", n3t);
                if (completeOutput) {
                    printTransition(fileOut, sif(), elemAux, a, ratestr);
                }
                if (!sit.goToElement(elemAux))
                    fileOut << "\n Error in transition";
                fileOut.flush();
                Proc.addRate(sif, sit, a, n3t * muCon);
            }

            elemAux.setCoordValue(EV, A1);
            sprintf(ratestr, "lamb%d", 1);
            if (completeOutput) {
                printTransition(fileOut, sif(), elemAux, a, ratestr);
            }
            if (!sit.goToElement(elemAux))
                fileOut << "\n Error in transition";
            fileOut.flush();
            Proc.addRate(sif, sit, a, lambClass1);

            elemAux.setCoordValue(EV, A2);
            sprintf(ratestr, "lamb%d", 2);
            if (completeOutput) {
                printTransition(fileOut, sif(), elemAux, a, ratestr);
            }
            if (!sit.goToElement(elemAux))
                fileOut << "\n Error in transition";
            fileOut.flush();
            Proc.addRate(sif, sit, a, lambClass2);

            elemAux.setCoordValue(EV, A3);
            sprintf(ratestr, "lamb%d", 3);
            if (completeOutput) {
                printTransition(fileOut, sif(), elemAux, a, ratestr);
            }
            if (!sit.goToElement(elemAux))
                fileOut << "\n Error in transition";
            fileOut.flush();
            Proc.addRate(sif, sit, a, lambClass3);
        }
    }

    cout << "\n Generating rewards... \n";
    if (completeOutput) {
        fileOut << "\n\nReward for State-actions:\n\n";
    }
    FLOAT auxValue, reward;
    for (sif.reset(); !sif.isDone(); sif.next()) {
        //n1 = sif().getCoordValue(N1);
        //n2 = sif().getCoordValue(N2);
        //n3 = sif().getCoordValue(N3);
        ev = sif().getCoordValue(EV);
        for (a = 0; a < Proc.getNActions(sif); ++a) {
            reactToAction(sif(), a, elemAux);
            n1t = elemAux.getCoordValue(N1);
            n2t = elemAux.getCoordValue(N2);
            n3t = elemAux.getCoordValue(N3);
            if (Proc.getStateRate(sif, a) == 0.0) {
                cout << "\n Error: zero State-action Rate:  "
                        << "state: [" << sif.getPosition() << "] " << sif()
                        << " action: " << a;
                cin >> ch;
                exit(1);
            }
            if (objFuncType == 1 || objFuncType == 3) {
                reward = (n1t * reward1 + n2t * reward2 + n3t * reward3) / Proc.getStateRate(sif, a);
                auxValue = UBReward / Proc.getStateRate(sif, a) - reward;
            } else {
                reward = (n1t * reward1 + n2t * reward2 + n3t * reward3) * muCon / Proc.getStateRate(sif, a);
                auxValue = UBReward * muCon / Proc.getStateRate(sif, a) - reward;
            }

            Proc.setCost(sif, a, auxValue);
            if (completeOutput) {
                fileOut << setw(2) << i << "   ";
                printState(fileOut, sif());
                fileOut << "   ";
                printAction(fileOut, sif(), a);
                fileOut << "   " << reward << "\n";
            }
        }
    }

    if (completeOutput) Proc.printOn(fileOut);
    else {
        fileOut << "\n\n*** Size of the constructed Continuos Time Markov Decision Process ***"
                << "\n\n  Number of states               : " << Proc.getNStates()
                << "\n  Number of state-actions        : " << Proc.getNStateActions()
                << "\n  Number of non-zero transitions : " << Proc.getNNoZerosTransitions()
                << "\n";
    }
    fileOut.flush();
    switch (optimizationMethod) {
        case 1:
        {
            CTMDPValueIterator values(Proc, policy, viPrecision,
                    tauFactor, maxIter);
            values.setNotModified();
            //values.setModified();
            values.setPrecision(viPrecision);
            if (useLastIterationData)
                values.setInitialValues(Val);
            values.go(maxIter);
            values.getValues(Val);
            fileOut << "\nClassical Values Iteration Algorithm"
                    << "\n   Optimal reward:       " << UBReward - values.getOptimalCost()
                    << "\n   Optimal reward error: " << values.getOptimalCostError()
                    << "\n   Number of iterations: " << values.getNumberOfIterations();

            //values.printResultsOn(fileOut);
        };
            break;
        case 2:
        {
            setFixedReservationPolicy(CS, policy, nWL,
                    nReservedWL1, nReservedWL2, nReservedWL3);
            fileOut << "\n\nGiven number of wavelenght for each class in node:";
            fileOut << "\n   Class 1 (quiting the ring in the node):  " << nReservedWL1;
            fileOut << "\n   Class 2 (crossing the node):             " << nReservedWL2;
            fileOut << "\n   Class 3 (entering the ring in the node): " << nReservedWL3;
            fileOut << "\n";
            for (sif.reset(); !sif.isDone(); sif.next()) {
                fileOut << "\n";
                printState(fileOut, sif());
                fileOut << "     ";
                printAction(fileOut, sif(), policy.getAction(sif));
                fileOut.flush();
            }
        };
            break;
        case 3:
        {
            CTVIANU values(Proc, policy, viPrecision, maxIter);
            values.setPrecision(viPrecision);
            values.setLogOn();
            values.setVerboseOn();
            values.go(maxIter);
            fileOut << "\nValues Iteration Algorithm without uniformization"
                    << "\n   Optimal reward:       " << UBReward - values.getOptimalCost()
                    << "\n   Optimal reward error: " << values.getOptimalCostError()
                    << "\n   Number of iterations: " << values.getNumberOfIterations();

            //values.printResultsOn(fileOut);
        };
            break;
        case 4:
        {
            CTVIAU values(Proc, policy, tauFactor, viPrecision, maxIter);
            values.setPrecision(viPrecision);
            values.setLogOn();
            values.setVerboseOn();
            values.go(maxIter);
            fileOut << "\nValues Iteration Algorithm without uniformization"
                    << "\n   Optimal reward:       " << UBReward - values.getOptimalCost()
                    << "\n   Optimal reward error: " << values.getOptimalCostError()
                    << "\n   Number of iterations: " << values.getNumberOfIterations();

            //values.printResultsOn(fileOut);
        };
            break;
        case 5:
        {
            CTVIASU values(Proc, policy, tauFactor, unifFactor,
                    viPrecision, maxIter);
            values.setPrecision(viPrecision);
            values.setLogOn();
            values.setVerboseOn();
            values.go(maxIter);
            fileOut << "\nValues Iteration Algorithm with semi-uniformization"
                    << "\n   Optimal reward:       " << UBReward - values.getOptimalCost()
                    << "\n   Optimal reward error: " << values.getOptimalCostError()
                    << "\n   Number of iterations: " << values.getNumberOfIterations();

            //values.printResultsOn(fileOut);
        };
            break;
        default:
        {
            fileOut << "\nInvalid algorithm type.";
            cout << "\nInvalid algorithm type.";
        };
    };
    CTMCEmbCTMDP chain(Proc, policy);
    CTMCStateProb stProb(chain);
    cout << " \n\n calculating limit probabilities...\n";
    //   if(useLastIterationData) {
    //	   for(sif.reset(); !sif.isDone(); sif.next()) {
    //         limitProb(sif.getPosition()) =
    //                  0.8*limitProb(sif.getPosition()) + 0.2/chain.getNStates();
    //      }
    //      stProb.setInitialProb(limitProb);
    //   }
    stProb.goToLimit(SWOR, lpPrecision, maxIter, limProbRelaxFactor);
    stProb.getStateProb(limitProb);



    //   Array<FLOAT> limitProb(Proc.getNStates(), 1.0/Proc.getNStates());
    //   for (sif.reset(); !sif.isDone(); sif.next())
    //      limitProb[sif.getPosition()] = stProb.getProb(sif);


    cout << "\n\n classifing states... \n";
    CTMkvChainStateClassifier classif(chain);
    classif.classify();
    cout << "\nFound:\n";
    if (classif.getNTransientStates() > 0) {
        cout << "\n   " << classif.getNTransientStates()
                << " transient state(s):";
    }
    cout << "\n   " << classif.getNClosedSets()
            << " closed irreductible state set(s):\n";

    if (completeOutput) {
        fileOut << "\n\nSTATIONARY POLICY AND LIMIT PROBABILITIES: \n";
    }
    FLOAT meanTimeBetweenDecisions = 0.0, activRate = 0.0,
            meanCustomers = 0.0, rejectionRate = 0.0;
    Array<FLOAT> NC1Prob(nWL + 1, 0.0), NC2Prob(nWL + 1, 0.0), NC3Prob(nWL + 1, 0.0);
    FLOAT throughput1 = 0.0, throughput2 = 0.0, throughput3 = 0.0,
            meanC1 = 0.0, meanC2 = 0.0, meanC3 = 0.0,
            meanReward1 = 0.0, meanReward2 = 0.0, meanReward3 = 0.0;

    const int TRANSORREC = 10000, MULTIREC = 20000;
    // NOCLASSIF and TRANSIENT are defined in ctclass.h

    bool ok;
    Array2<int> SClassif1(nWL + 1, nWL + 1, NOCLASSIF),
            SClassif2(nWL + 1, nWL + 1, NOCLASSIF),
            SClassif3(nWL + 1, nWL + 1, NOCLASSIF);
    int sC, sCRec, sifC, newSC;
    for (sif.reset(); !sif.isDone(); sif.next()) {
        n1 = sif().getCoordValue(N1);
        n2 = sif().getCoordValue(N2);
        n3 = sif().getCoordValue(N3);
        ev = sif().getCoordValue(EV);
        a = policy.getAction(sif);
        auxValue = limitProb[sif.getPosition()];
        reactToAction(sif(), a, elemAux);
        n1t = elemAux.getCoordValue(N1);
        n2t = elemAux.getCoordValue(N2);
        n3t = elemAux.getCoordValue(N3);
        NC1Prob(n1t) += auxValue;
        NC2Prob(n2t) += auxValue;
        NC3Prob(n3t) += auxValue;
        meanC1 += n1t*auxValue;
        meanC2 += n2t*auxValue;
        meanC3 += n3t*auxValue;
        meanReward1 += reward1 * n1t*auxValue;
        meanReward2 += reward2 * n2t*auxValue;
        meanReward3 += reward3 * n3t*auxValue;
        throughput1 += n1t * muCon*auxValue;
        throughput2 += n2t * muCon*auxValue;
        throughput3 += n3t * muCon*auxValue;
        sifC = classif.getStateClassification(sif);
        if (completeOutput) {
            printState(fileOut, sif());
            fileOut << "     ";
            printAction(fileOut, sif(), a);
            fileOut << "    ";
            if (sifC == TRANSIENT)
                fileOut << "T " << "     ";
            else
                fileOut << "R" << sifC << "     ";
            fileOut << auxValue << "\n";
            cout.flush();
        }

        if (ev == A1 || ev == A2 || ev == A3) {
            switch (ev) {
                case A1:
                    sC = SClassif1(n2, n3);
                    break;
                case A2:
                    sC = SClassif2(n1, n3);
                    break;
                case A3:
                    sC = SClassif3(n1, n2);
                    break;
            }
            if (sC != sifC && sC != MULTIREC) {
                ok = false;
                if (sC == NOCLASSIF) {
                    newSC = sifC;
                    ok = true;
                }
                if (!ok && sC == TRANSIENT) {
                    newSC = TRANSORREC + sifC;
                    ok = true;
                }
                if (!ok && sifC == TRANSIENT) {
                    if (sC > TRANSORREC)
                        newSC = sC;
                    else
                        newSC = TRANSORREC + sC;
                    ok = true;
                }
                if (!ok && sC < TRANSORREC) {
                    if (sC != sifC)
                        newSC = MULTIREC;
                    ok = true;
                }
                if (!ok && sC > TRANSORREC) {
                    if (sC != TRANSORREC + sifC)
                        newSC = MULTIREC;
                    ok = true;
                }
                if (ok)
                    switch (ev) {
                        case A1:
                            SClassif1(n2, n3) = newSC;
                            break;
                        case A2:
                            SClassif2(n1, n3) = newSC;
                            break;
                        case A3:
                            SClassif3(n1, n2) = newSC;
                            break;
                    }
            }
        }
    }

    fileOut << "\n\n\nState classification in Class 1 arrival periods:";
    fileOut << "\n\nn2\\n3";
    for (j = 0; j <= nWL; ++j)
        fileOut << "  " << setw(2) << j;
    fileOut << "\n";
    for (i = 0; i <= nWL; ++i) {
        fileOut << "\n" << setw(2) << i << ":  ";
        for (j = 0; j <= nWL; ++j) {
            fileOut << "  ";
            sC = SClassif1(i, j);
            ok = false;
            //         fileOut << setw(2) << sC;
            switch (sC) {
                case NOCLASSIF:
                    fileOut << " .";
                    ok = true;
                    break;
                case TRANSIENT:
                    fileOut << " T";
                    ok = true;
                    break;
                case MULTIREC:
                    fileOut << "mr";
                    ok = true;
                    break;
            }
            if (!ok && sC < TRANSORREC) {
                fileOut << "R" << setw(1) << sC;
                ok = true;
            }
            if (!ok && sC < MULTIREC) {
                fileOut << "r" << setw(1) << sC - TRANSORREC;
                ok = true;
            }
            if (!ok)
                fileOut << "??";
        }
    }
    fileOut << "\n\n\nState classification in Class 2 arrival periods:";
    fileOut << "\n\nn1\\n3";
    for (j = 0; j <= nWL; ++j)
        fileOut << "  " << setw(2) << j;
    fileOut << "\n";
    for (i = 0; i <= nWL; ++i) {
        fileOut << "\n" << setw(2) << i << ":  ";
        for (j = 0; j <= nWL; ++j) {
            fileOut << "  ";
            sC = SClassif2(i, j);
            ok = false;
            //         fileOut << setw(2) << sC;
            switch (sC) {
                case NOCLASSIF:
                    fileOut << " .";
                    ok = true;
                    break;
                case TRANSIENT:
                    fileOut << " T";
                    ok = true;
                    break;
                case MULTIREC:
                    fileOut << "mr";
                    ok = true;
                    break;
            }
            if (!ok && sC < TRANSORREC) {
                fileOut << "R" << setw(1) << sC;
                ok = true;
            }
            if (!ok && sC < MULTIREC) {
                fileOut << "r" << setw(1) << sC - TRANSORREC;
                ok = true;
            }
            if (!ok)
                fileOut << "??";
        }
    }
    fileOut << "\n\n\nState classification in Class 3 arrival periods:";
    fileOut << "\n\nn1\\n2";
    for (j = 0; j <= nWL; ++j)
        fileOut << "  " << setw(2) << j;
    fileOut << "\n";
    for (i = 0; i <= nWL; ++i) {
        fileOut << "\n" << setw(2) << i << ":  ";
        for (j = 0; j <= nWL; ++j) {
            fileOut << "  ";
            sC = SClassif3(i, j);
            ok = false;
            //         fileOut << setw(2) << sC;
            switch (sC) {
                case NOCLASSIF:
                    fileOut << " .";
                    ok = true;
                    break;
                case TRANSIENT:
                    fileOut << " T";
                    ok = true;
                    break;
                case MULTIREC:
                    fileOut << "mr";
                    ok = true;
                    break;
            }
            if (!ok && sC < TRANSORREC) {
                fileOut << "R" << setw(1) << sC;
                ok = true;
            }
            if (!ok && sC < MULTIREC) {
                fileOut << "r" << setw(1) << sC - TRANSORREC;
                ok = true;
            }
            if (!ok)
                fileOut << "??";
        }
    }

    /*
            fileOut << "\n\nSTATIONARY POLICY RESERVES WAKELENGHT IN THE STATES: \n";
       int max, min;
            for(sif.reset(); !sif.isDone(); sif.next()) {
          n1 = sif().getCoordValue(N1);
          n2 = sif().getCoordValue(N2);
          n3 = sif().getCoordValue(N3);
          ev = sif().getCoordValue(EV);
          if(Proc.getNActions(sif)==2 && policy.getAction(sif)==REJ) {
             printState(fileOut, sif());
               fileOut << "     ";
             printAction(fileOut, sif(), a);
             fileOut << "    ";
             if(classif.getStateClassification(sif)==TRANSIENT)
                fileOut << "T";
             else
                fileOut << "R";
             max = ((n1+n2)>=(n2+n3))? (n1+n2) : (n2+n3);
             min = ((n1+n2)<=(n2+n3))? (n1+n2) : (n2+n3);
             fileOut << "      n1+n2 = " << setw(2) << (n1+n2)
                        << "   n2+n3 = " << setw(2) << (n2+n3)
                        << "   max = " << setw(2) << max
                        << "   min = " << setw(2) << min;
             fileOut << "\n";
                    cout.flush();
          }
            }

            Policy minNWLPolicy(CS);
       INDEX minNWL1, minNWL2, minNWL3;
       calcMinReservationLevelsForAPolicy(CS, policy, nWL,
                                                     minNWL1, minNWL2, minNWL3);
            setFixedReservationPolicy(CS, minNWLPolicy, nWL,
                                                     minNWL1, minNWL2, minNWL3);
       bool isDiff=false;
            fileOut << "\n\nTHE POLICY IS THE RESERVED number of wavelenght for each class: \n";
       fileOut << "\n   Class 1 (quiting the ring in the node):  " << minNWL1;
       fileOut << "\n   Class 2 (crossing the node):             " << minNWL2;
       fileOut << "\n   Class 3 (entering the ring in the node): " << minNWL3;
            for(sif.reset(); !sif.isDone(); sif.next()) {
          if(policy.getAction(sif)!=minNWLPolicy.getAction(sif)) {
             if(!isDiff) {
                    fileOut << "\n\nEXCEPT FOR THE STATES: \n";
                isDiff = true;
             }
             printState(fileOut, sif());
               fileOut << "     ";
             printAction(fileOut, sif(), a);
             fileOut << "    ";
             if(classif.getStateClassification(sif)==TRANSIENT)
                fileOut << "T";
             else
                fileOut << "R";
             fileOut << "\n";
                    cout.flush();
          }
            }
     */
    PolicyAnalysis(fileOut, CS, policy, nWL);

    //fileOut << "\n\nPERIODO PARA CADA CONJUNTO FECHADO DE ESTADOS:";
    //for(int cs=0; cs<classif.getNClosedSets(); ++cs)
    //  fileOut << "\nset R" << (cs+1) << ":  "
    //          << classif.getClosedSetPeriod(cs);

    fileOut << "\n\nLimit Probability for number of clients in each Class:";
    for (n = 0; n <= nWL; ++n)
        fileOut << "\n   " << n << ": " << setw(15) << NC1Prob(n)
        << "     " << setw(15) << NC2Prob(n)
        << "     " << setw(15) << NC3Prob(n);
    fileOut << "\n\nMean number of calls of each class:";
    fileOut << "\n   Class 1: " << setw(15) << meanC1
            << "\n   Class 2: " << setw(15) << meanC2
            << "\n   Class 3: " << setw(15) << meanC3;

    fileOut << "\n\nThroughput of each class:";
    fileOut << "\n   Class 1: " << setw(15) << throughput1
            << "\n   Class 2: " << setw(15) << throughput2
            << "\n   Class 3: " << setw(15) << throughput3;
    probBlock1 = 1.0 - throughput1 / lambClass1;
    probBlock2 = 1.0 - throughput2 / lambClass2;
    probBlock3 = 1.0 - throughput3 / lambClass3;
    fileOut << "\n\nBlock probability of each class:";
    fileOut << "\n   Class 1: " << setw(15) << probBlock1
            << "\n   Class 2: " << setw(15) << probBlock2
            << "\n   Class 3: " << setw(15) << probBlock3;

    fileOut << "\n\nMean reward of each class:";
    fileOut << "\n   Class 1: " << setw(15) << meanReward1
            << "\n   Class 2: " << setw(15) << meanReward2
            << "\n   Class 3: " << setw(15) << meanReward3
            << "\n\n   Total:   " << setw(15) << meanReward1 +
            meanReward2 + meanReward3;
}

void Store_Matrix_lambCon(int numNos) {
    for (int i = 0; i < numNos; i++) {
        for (int j = 0; j < numNos; j++) {
            Matrix_lambCon[i][j] = lambCon;
        }//fim for
    }//fim for
}//fim m�todo

int indice(int i, int j) {
    if (i + j < 0) {
        return i + j + numNos;
    } else if (i + j > numNos) {
        return i + j - numNos;
    } else
        return i + j;
}

void Set_lambClass(int noRef, int numNos) {
    //Taxas de cada classe para o primeiro n� e considerando anel de 4 n�s unidirecional;

    lambClass1 = Matrix_lambCon[indice(noRef, -1)][noRef]*(1 - Set_ProbBlock[indice(noRef, -1)][indice(noRef, -1)][noRef]) +
            Matrix_lambCon[indice(noRef, -2)][noRef]*((1 - Set_ProbBlock[indice(noRef, -2)][indice(noRef, -2)][noRef])*(1 - Set_ProbBlock[indice(noRef, -1)][indice(noRef, -2)][noRef])) +
            Matrix_lambCon[indice(noRef, -3)][noRef]*((1 - Set_ProbBlock[indice(noRef, -3)][indice(noRef, -3)][noRef])*(1 - Set_ProbBlock[indice(noRef, -2)][indice(noRef, -3)][noRef])*(1 - Set_ProbBlock[indice(noRef, -1)][indice(noRef, -3)][noRef]));
    cout << "lambClass1: " << lambClass1 << "\n";

    lambClass2 = Matrix_lambCon[indice(noRef, -1)][indice(noRef, 1)]*(1 - Set_ProbBlock[indice(noRef, -1)][indice(noRef, -1)][indice(noRef, 1)]) +
            Matrix_lambCon[indice(noRef, -1)][indice(noRef, 2)]*((1 - Set_ProbBlock[indice(noRef, -1)][indice(noRef, -1)][indice(noRef, 2)])*(1 - Set_ProbBlock[indice(noRef, 1)][indice(noRef, -1)][indice(noRef, 2)])) +
            Matrix_lambCon[indice(noRef, -2)][indice(noRef, 1)]*((1 - Set_ProbBlock[indice(noRef, -2)][indice(noRef, -2)][indice(noRef, 1)])*(1 - Set_ProbBlock[indice(noRef, -1)][indice(noRef, -2)][indice(noRef, 1)]));
    cout << "lambClass2: " << lambClass2 << "\n";

    lambClass3 = Matrix_lambCon[noRef][indice(noRef, 1)] + Matrix_lambCon[noRef][indice(noRef, 2)] + Matrix_lambCon[noRef][indice(noRef, 3)];
    cout << "lambClass3: " << lambClass3 << "\n";
    //system("pause");

}//fim m�todo

class Node {
public:
    void Box(int);
    //void Matrix_ProbBlock(int);
    void Update_Set_ProbBlock(int, float);
};

void Node::Update_Set_ProbBlock(int noRef, float probBlock1, float probBlock2, float probBlock3) {
    if (noRef == 1) {
        Set_ProbBlock[noRef][i][j] = probBlock1;
        ...//Implementar as coisas aqui..
    } else if (noRef == 2) {
        Set_ProbBlock[noRef][i][j] = probBlock3;
        ...
    } else if (noRef == 3) {
        Set_ProbBlock[noRef][i][j] = probBlock2;
        ...
    } else if (noRef == 4) {
        Set_ProbBlock[noRef][i][j] = probBlock3;
        ...
    }
}

void Node::Box(int numNos) {
    for (int k = 0; k < numNos; k++) // varre os frames
    {
        for (int i = 0; i < numNos; i++) {
            for (int j = 0; j < numNos; j++) {

                Set_ProbBlock[k][i][j] = 0;
            }//fim for
        }//fim for

    }//fim for
    //        system("pause");
}//fim m�todo

void Node::Update_ProbBlock(int noRef) {

}//fim m�todo

void main(int argc, char *argv[]) {
    IDTimer crono;
    char ch, name[20], nameIn[100], nameOut[100];
    cout << "Input and output file name: ";
    cin >> name;
    cout << "\n";
    strcpy(nameIn, PATH);
    strcat(nameIn, name);
    strcat(nameIn, ".in");
    getDataFromFile(nameIn);
    strcpy(nameOut, PATH);
    strcat(nameOut, name);
    strcat(nameOut, ".out");
    ofstream fileOut(nameOut);

    showDataOn(fileOut);

    MultiLinearSet S(4);
    S.setMinMaxCoordValue(N1, 0, nWL);
    S.setMinMaxCoordValue(N2, 0, nWL);
    S.setMinMaxCoordValue(N3, 0, nWL);
    S.setMinMaxCoordValue(EV, 0, 5);

    OpticalWDMCondition cond(nWL);


    ConstrMultiLinearSet CS(S, cond);

    FLOAT //probLoss1=0.0, probLoss2=0.0, probLoss3=0.0,
    probBlock1, probBlock2, probBlock3;
    //  FLOAT newProbLoss1,  newProbLoss2,  newProbLoss3;
    FLOAT IncomingNetRateInNode, OutgoingNetRateInNode,
            CrossingInNodeRate, CrossingOutNodeRate;
    FLOAT error = 1.0, aux, iterTime, totalTime = 0.0;
    Array2<FLOAT> lambIn(4, 4, 0.0), lambOut(4, 4, 0.0),
            lambInInNextNode(4, 4, 0.0);
    if (useLastIterationData) {
        Val.restart(CS.getNElements(), 0.0);
        limitProb.restart(CS.getNElements(), 1.0 / CS.getNElements());
    }

    Node no; //primeira inst�ncia do objeto da classe node
    int k = 0; // primeiro frame da matriz
    int noRef = 0;

    no.Box(numNos); //inicializa o box -> conjunto de frames
    //no.Matrix_ProbBlock(k); //

    Store_Matrix_lambCon(numNos);

    no.Update_Set_ProbBlock(noRef, probBlock1, probBlock2, probBlock3);
    //tem q ver onde chamar este metodo
    /*
      lambIn(0, 1) = lambCon;
      lambIn(0, 2) = lambCon;
      lambIn(0, 3) = lambCon;

      lambIn(1, 0) = lambCon;
      lambIn(2, 0) = lambCon;
      lambIn(3, 0) = lambCon;

      lambIn(2, 1) = lambCon;
      lambIn(3, 1) = lambCon;
      lambIn(3, 2) = lambCon;
     */

    int iteration = 0;
    while (error > lpPrecision) { // (error>0.1*CS.getNElements()*lpPrecision) {
        crono.start();
        iterTime = 0.0;
        Set_lambClass(noRef, numNos);
        /*
             lambClass1 = lambIn(1, 0) + lambIn(2, 0) + lambIn(3, 0);
             lambClass2 = lambIn(2, 1) + lambIn(3, 1) + lambIn(3, 2);
             lambClass3 = lambIn(0, 1) + lambIn(0, 2) + lambIn(0, 3);
         */
        ++iteration;
        fileOut << "\n\n************ ITERATION " << iteration << " ************";
        cout << "\n\n************ ITERATION " << iteration << " ************\n\n";
        if (objFuncType == 3 || objFuncType == 4) {
            cout << "\nClass 1 reward: " << reward1;
            cout << "\nClass 2 reward: " << reward2;
            cout << "\nClass 3 reward: " << reward3;
            fileOut << "\nClass 1 reward: " << reward1;
            fileOut << "\nClass 2 reward: " << reward2;
            fileOut << "\nClass 3 reward: " << reward3;
        }
        fileOut << "\n";

        go(CS, fileOut, lambClass1, lambClass2, lambClass3,
                probBlock1, probBlock2, probBlock3);
        /*
            lambOut(0, 1) = (1.0-probBlock3)*lambIn(0, 1);
            lambOut(0, 2) = (1.0-probBlock3)*lambIn(0, 2);
            lambOut(0, 3) = (1.0-probBlock3)*lambIn(0, 3);

            lambOut(1, 0) = (1.0-probBlock1)*lambIn(1, 0);
            lambOut(2, 0) = (1.0-probBlock1)*lambIn(2, 0);
            lambOut(3, 0) = (1.0-probBlock1)*lambIn(3, 0);

            lambOut(2, 1) = (1.0-probBlock2)*lambIn(2, 1);
            lambOut(3, 1) = (1.0-probBlock2)*lambIn(3, 1);
            lambOut(3, 2) = (1.0-probBlock2)*lambIn(3, 2);

            IncomingNetRateInNode = lambIn(0, 1)  + lambIn(0, 2)  + lambIn(0, 3);
            OutgoingNetRateInNode = lambOut(1, 0) + lambOut(2, 0) + lambOut(3, 0);
            CrossingInNodeRate    = lambIn(2, 1)  + lambIn(3, 1)  + lambIn(3, 2);
            CrossingOutNodeRate   = lambOut(3, 2) + lambOut(0, 2) + lambOut(0, 3);

                 fileOut << "\n\nRESULTS FOR HOMOGENEOUS 4 NODE NETWORK:";
            fileOut  << "\n\nprobability of blocking a class 1 arriving user: " << probBlock1;
            fileOut  << "  \nprobability of blocking a class 2 arriving user: " << probBlock2;
            fileOut  << "  \nprobability of blocking a class 3 arriving user: " << probBlock3;
                 fileOut << "\n\nprobability of blocking";
                 fileOut << "\n   a 2 nodes call: " << 1.0-( (1.0-probBlock3)*
                                                        (1.0-probBlock1)  )
                         << "\n   a 3 nodes call: " << 1.0-( (1.0-probBlock3)*
                                                        (1.0-probBlock2)*
                                                        (1.0-probBlock1)  )
                         << "\n   a 4 nodes call: " << 1.0-( (1.0-probBlock3)*
                                                        (1.0-probBlock2)*
                                                        (1.0-probBlock2)*
                                                        (1.0-probBlock1)  );
                 fileOut << "\n\nArriving in the Network by Node 1 Rate: " << IncomingNetRateInNode;
                 fileOut << "  \nGoing   out the Network by Node 1 Rate: " << OutgoingNetRateInNode;
                 fileOut << "  \nCrossing the Node 1 Rate (in):      " << CrossingInNodeRate;
                 fileOut << "  \nCrossing the Node 1 Rate (out):     " << CrossingOutNodeRate;

            fileOut << "\n\n";
                 fileOut << "\n\nArrival rate in node 1:";
            for(int i=0; i<4; ++i) {
                    fileOut << "\n" << (i+1) << ":";
               for(int j=0; j<4; ++j) {
                  fileOut  << "  " << setw(15) << lambIn(i, j);
               }
            }

                 fileOut << "\n\nOutput rate in node 1:";
            for(int i=0; i<4; ++i) {
                    fileOut << "\n" << (i+1) << ":";
               for(int j=0; j<4; ++j) {
                   fileOut  << "  " << setw(15) << lambOut(i, j);
               }
            }

                 fileOut << "\n\nBlock Probability in node 1:";
            for(int i=0; i<4; ++i) {
                    fileOut << "\n" << (i+1) << ":";
               for(int j=0; j<4; ++j) {
                  if(lambIn(i, j)!=0.0)
                     fileOut  << "  " << setw(15)
                              << (lambIn(i, j)-lambOut(i, j))/lambIn(i, j);
                  else
                     fileOut  << "  ---------------";
               }
            }

               lambInNextNodeClass1 = lambOut(2, 1) + lambOut(3, 1) + lambOut(0, 1);
               lambInNextNodeClass2 = lambOut(3, 2) + lambOut(0, 2) + lambOut(0, 3);
               lambInNextNodeClass3 = 3*lambCon;


                 fileOut << "\n\nnew arrival rate in the next node:";
                 fileOut << "\n   Class 1: " << setw(15) << lambInNextNodeClass1
                         << "\n   Class 2: " << setw(15) << lambInNextNodeClass2
                         << "\n   Class 3: " << setw(15) << lambInNextNodeClass3;
               fileOut << "\n\n";
         */
        error = fabs(lambInNextNodeClass1 - lambClass1);
        aux = fabs(lambInNextNodeClass2 - lambClass2);
        if (aux > error) error = aux;
        aux = fabs(lambInNextNodeClass3 - lambClass3);
        if (aux > error) error = aux;
        fileOut << "\n\nERROR IN ITERATION " << iteration
                << ": " << setw(15) << error << "\n\n";
        cerr << "\n\nERROR IN ITERATION " << iteration
                << ": " << setw(15) << error << "\n\n";
        /*
        // preparing for the next iteration:

        lambIn(0, 1) = lambCon;
        lambIn(0, 2) = lambCon;
        lambIn(0, 3) = lambCon;

        lambIn(1, 0) = (1.0-rateUpdateRelaxFactor)*lambIn(1, 0) +
                                     rateUpdateRelaxFactor*lambOut(2, 1);
        lambIn(2, 0) = (1.0-rateUpdateRelaxFactor)*lambIn(2, 0) +
                                     rateUpdateRelaxFactor*lambOut(3, 1);
        lambIn(3, 0) = (1.0-rateUpdateRelaxFactor)*lambIn(3, 0) +
                                     rateUpdateRelaxFactor*lambOut(0, 1);

        lambIn(2, 1) = (1.0-rateUpdateRelaxFactor)*lambIn(2, 1) +
                                     rateUpdateRelaxFactor*lambOut(3, 2);
        lambIn(3, 1) = (1.0-rateUpdateRelaxFactor)*lambIn(3, 1) +
                                     rateUpdateRelaxFactor*lambOut(0, 2);
        lambIn(3, 2) = (1.0-rateUpdateRelaxFactor)*lambIn(3, 2) +
                                     rateUpdateRelaxFactor*lambOut(0, 3);

        if (objFuncType==3 || objFuncType==4) {
/*
           reward1 = ((1.0-probBlock3)*
                      (1.0-probBlock2)*
                      (1.0-probBlock2)*
                      (1.0           )*lambOut(2, 1) +
                      (1.0-probBlock3)*
                      (1.0-probBlock2)*
                      (1.0           )*lambOut(3, 1) +
                      (1.0-probBlock3)*
                      (1.0            )*lambOut(0, 1)  )/(lambOut(2, 1) +
                                                         lambOut(3, 1) +
                                                         lambOut(0, 1)      );
           reward2 = ((1.0-probBlock3)*
                      (1.0-probBlock2)*
                      (1.0           )*
                      (1.0-probBlock1)*lambOut(2, 1) +
                      (1.0-probBlock3)*
                      (1.0           )*
                      (1.0-probBlock1)*lambOut(3, 1) +
                      (1.0-probBlock3)*
                      (1.0           )*
                      (1.0-probBlock2)*
                      (1.0-probBlock1)*lambOut(3, 2)   )/(lambOut(2, 1) +
                                                          lambOut(3, 1) +
                                                          lambOut(3, 2)    );
           reward3 = ((1.0           )*
                      (1.0-probBlock1)*lambOut(0, 1) +
                      (1.0           )*
                      (1.0-probBlock2)*
                      (1.0-probBlock1)*lambOut(0, 2) +
                      (1.0           )*
                      (1.0-probBlock2)*
                      (1.0-probBlock2)*
                      (1.0-probBlock1)*lambOut(0, 3)   )/(lambOut(0, 1) +
                                                          lambOut(0, 2) +
                                                          lambOut(0, 3)    );
           reward1 = ( lambOut(3, 0)/2 +
                       lambOut(2, 0)/3 +
                       lambOut(1, 0)/4       )/(lambOut(3, 0) +
                                              lambOut(2, 0) +
                                              lambOut(1, 0)      );
           reward2 = ((1.0-probBlock1)*lambOut(2, 1)/4 +
                      (1.0-probBlock1)*lambOut(3, 1)/3 +
                      (1.0-probBlock2)*
                      (1.0-probBlock1)*lambOut(3, 2)/4   )/(lambOut(2, 1) +
                                                          lambOut(3, 1) +
                                                          lambOut(3, 2)    );
           reward3 = ((1.0-probBlock1)*lambOut(0, 1)/2 +
                      (1.0-probBlock2)*
                      (1.0-probBlock1)*lambOut(0, 2)/3 +
                      (1.0-probBlock2)*
                      (1.0-probBlock2)*
                      (1.0-probBlock1)*lambOut(0, 3)/4   )/(lambOut(0, 1) +
                                                          lambOut(0, 2) +
                                                          lambOut(0, 3)    );

         */ /*
           reward1 = ( lambIn(3, 0)/2 +
                       lambIn(2, 0)/3 +
                       lambIn(1, 0)/4       )/( lambIn(3, 0) +
                                                lambIn(2, 0) +
                                                lambIn(1, 0)      );
           reward2 = ((1.0-probBlock1)* lambIn(2, 1)/4 +
                      (1.0-probBlock1)* lambIn(3, 1)/3 +
                      (1.0-probBlock2)*
                      (1.0-probBlock1)* lambIn(3, 2)/4   )/( lambIn(2, 1) +
                                                             lambIn(3, 1) +
                                                             lambIn(3, 2)    );
           reward3 = ((1.0-probBlock1)* lambIn(0, 1)/2 +
                      (1.0-probBlock2)*
                      (1.0-probBlock1)* lambIn(0, 2)/3 +
                      (1.0-probBlock2)*
                      (1.0-probBlock2)*
                      (1.0-probBlock1)* lambIn(0, 3)/4   )/( lambIn(0, 1) +
                                                             lambIn(0, 2) +
                                                             lambIn(0, 3)    );
           if (objFuncType==4) {
              reward1 = muCon*reward1;
              reward2 = muCon*reward2;
              reward3 = muCon*reward3;
           }
           //reward1 /= reward3;
           //reward2 /= reward3;
           //reward3 = 1.0;
        }*/
        crono.mark();
        iterTime = crono.elapsedTime();
        totalTime += iterTime;
        fileOut << "\n\nTime of iteration: " << iterTime << "\n\n";
        cout << "\n\nTime of iteration: " << iterTime << "\n\n";
        FLOAT netThroughput, JainFairness, soma, soma2, auxF;
        netThroughput = 4 * OutgoingNetRateInNode;
        auxF = (1.0 - probBlock3)*(1.0 - probBlock1);
        soma = 4 * auxF;
        soma2 = 4 * auxF*auxF;
        auxF = (1.0 - probBlock3)*(1.0 - probBlock2)*(1.0 - probBlock1);
        soma += 4 * auxF;
        soma2 += 4 * auxF*auxF;
        auxF = (1.0 - probBlock3)*(1.0 - probBlock2)*(1.0 - probBlock2)*(1.0 - probBlock1);
        soma += 4 * auxF;
        soma2 += 4 * auxF*auxF;
        JainFairness = soma * soma / (12 * soma2);
        cerr << "\n\nThroughput of the 4-nodes network:       " << netThroughput;
        cerr << "\nJain's Fairness for the 4-nodes network: " << JainFairness << "\n\n";
        fileOut << "\n\nThroughput of the 4-nodes network:       " << netThroughput;
        fileOut << "\nJain's Fairness for the 4-nodes network: " << JainFairness << "\n\n";
    }



    fileOut << "\n\nTOTAL TIME: " << totalTime << "\n\n";
    cout << "\n\nTOTAL TIME: " << totalTime << "\n\n";

    fileOut.flush();
    cout << "\nThat's all...\n";
    cin >> ch;
}






